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We use a simple fragmentation model to describe the statistical behavior of the Voronoi cell 
patterns generated by a homogeneous and isotropic set of points in ID and in 2D. In particular, we 
are interested in the distribution of sizes of these Voronoi cells. Our model is completely defined 
by two probability distributions in ID and again in 2D, the probability to add a new point inside 
an existing cell and the probability that this new point is at a particular position relative to the 
preexisting point inside this cell. In ID the first distribution depends on a single parameter while the 
second distribution is defined through a fragmentation kernel; in 2D both distributions depend on 
a single parameter. The fragmentation kernel and the control parameters are closely related to the 
physical properties of the specific system under study. We use our model to describe the Voronoi cell 
patterns of several systems. Specifically, we study the island nucleation with irreversible attachment, 
the ID car-parking problem, the formation of second-level administrative divisions, and the pattern 
formed by the Paris Metro stations. 

PACS numbers: 89.75.Kd,68.55.A-,05.40.-a,81.16.Rf 



I. INTRODUCTION 

Consider a set of points — usually called centers [1-3] 
even though they are not geometric centers — on a dis- 
crete lattice. The Voronoi cell of a particular center i, is 
defined by all lattice points which are closer to i than any 
other center. Figure 1 shows a typical two-dimensional 
Voronoi pattern for the case where the positions of the 
centers are completely uncorr elated. This case is usually 
called Poisson Voronoi (PV). When the position of the 
centers are correlated, the system is called non-Poisson 
Voronoi (NPV). 

Many different systems in nature resemble the PV pat- 
terns. Some examples can be found in areas such as ecol- 
ogy, astronomy, geology, biology, physics, and meteorol- 
ogy; see Refs. [1-6]. Applications of NPV patterns are 
not so extensive as those of the PV case. A couple of ex- 
amples of NPV applications can be found in Refs. [3, 7-9]. 

One of the most important quantities in this system 
is the distribution of the sizes of the Voronoi cells P{S). 
The scaled size is defined a>s s = S/ (S), where (S) is 
the average of S. The scaled size distribution is given by 
P{s) = (S) P{s (S)). There are many theoretical and nu- 
merical studies about PV systems [1, 2, 10-15]. In spite of 
this, the patterns formed by Voronoi cells have not been 
understood completely even in this simplest case, where 
the positions of the centers are not correlated. Most of 
our knowledge is based on empirical equations and nu- 
merical simulations. In fact, an analytical expression for 
P{s) is known just for the ID case with uncorrelated 
centers [10, 11]. 

In this paper we focus on NPV patterns in ID and 2D. 
In Sec. II we review some important properties of the PV 
cells. In Sec. Ill we propose a model to generate NPV 
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FIG. 1. Typical pattern of Poisson Voronoi cells. Note that in 
general the positions of the centers inside cells do not coincide 
with the geometrical centers of cells. 

cells for a homogeneous and isotropic set of points. In 
Sec. IV we provide several examples of different systems 
which can be described by NPV patterns. Finally, in Sec. 
V we give conclusions. 



II. POISSON VORONOI CELLS 

A. One-dimensional Poisson Voronoi cells 

In one dimension, we have a ring divided in several 
sections called gaps. For the PV case, the positions 
of centers are completely random. Then, the probabil- 
ity density to find a gap with a length between X and 
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X + dX^ p^^\X)^ is given by pe~^^, where p is the 
density of centers. The normahzed gap size is defined 
as X = X/{X) where (X) = is the average gap size. 
Thus, the normahzed gap size distribution can be written 
as p^^\s) = . By definition, p^'^\x) is the probabil- 
ity density to find a gap which starts and ends with a 
center subject to the condition that there are n addi- 
tional centers inside the gap. These distributions sat- 
isfy the normalization conditions dxp^'^\x) = 1 and 
dxxp^^\x) = n + 1. Because sizes of the adjacent 
gaps are not correlated, it is possible to write the nor- 
malized spacing distributions for arbitrary values of n in 
Laplace space as 



p(")(Z) 



(p(°)(Z)) 



n+1 



(1) 



where p^^^ (/) = dxe ^ ^ p^'^'^ (x) is the Laplace trans- 
form of p^^\x) [16]. Consequently, the normalized spac- 
ing distributions between centers are given by 



The pair correlation function has the simple form 

oo 

5(x) = ^p(")(;r) = l. 
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The size distribution of the Vononoi cells is related to the 
next-nearest-neighbor distribution according to [8] 



P(s) = 2p(i)(2s). 



(4) 



This equation is a consequence of the one- dimensional na- 
ture of the ring. Therefore, the simple relation between 
P{s) and p^^\s) shown above is not valid for higher di- 
mensions. Explicitly, the distribution of sizes of Voronoi 
cells in ID is 
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The next radial distribution p^^\R) is given by 
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P^^\R) = 27r pR / dR'p^''\R')Q{R',R), (8) 
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where Q{R' , R) = e-^^^^ ) is the probability density 
to have the annulus R' < r < R free of centers. The 
integral in Eq. (8) can be calculated straightforwardly to 
give 



p^^\R)=2^^ p" R^e-^P""". 



(9) 



Following the previous procedure, one can show that, for 
arbitrary n, 

n\ 

As for multineighbor spacing distributions in ID [16- 
18], the radial distributions have information about the 
structure of the system. For example, since 
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it is possible to calculate the radial distribution function 
g{R). From Eqs. (10) and (11) it is easy to find 
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This result is consistent with this case of uncorrelated 
centers g{R) = 1. In general the concentration of centers 
c{R) can be extracted from p^^\R). From Eq. (6) it 
follows 
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B. Two-dimensional Poisson Voronoi cells 

Let p^^\R) be the radial probability density that, 
given an island at r = 0, its (n + 1)^^ neighbor is be- 
tween R and R + dR. Given that there are n additional 
centers inside of the circle of radius i?, the radial spacing 
distribution, p^^\R)^ can be calculated as follows. As 
usual, p is the density of centers. On average, the total 
number of centers, c(i^), within a disk with radius R is 
27r pR. It is weh known [19] that c{R) and p^^\R) are 
related by 



^(o)(i^) =c(i?)e- 



- dr 27rr c(r) 
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More details about this equation are given in the Ap- 
pendix. Since c{R) = 2 7rRp, p^^\R) has the simple 
form 
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Even though we know the exact expression for the ra- 
dial spacing distributions for arbitrary values of n in the 
PV case, the exact functional form of P (5) is not known 
for d > 2. In fact, just a few exact analytical results are 
known in this case. One of them was reported in 1962 
by Gilbert [15], who showed that the second moment of 
P{S) is 0.280(5'). 

As mentioned previously, obtaining the exact expres- 
sion of P (5) for d > 1 is quite complicated partially due 
geometrical complications. In d = 1 each new center di- 
vides just one of the existing Voronoi cells to form two 
new cells. This fact leads to the simple relation between 
P{x) and p^^\x) given by Eq. (4). In higher dimensions 
this does not apply; the Voronoi cell of a new center 
is formed at the expense of several preexisting Voronoi 
cells. An analogous relation to Eq. (4) ioi d> 1 remains 
unknown, and it could involve several p^^^ (r) in a non- 
trivial way. However, it is well accepted that P{s) can be 
approximated by the gamma distribution 110,(5). Based 
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FIG. 2. (Color online) Poisson Voronoi cell-size distribution 
P(s) with d — 2. The agreement between the numerical re- 
sults and Eq. (14) is excellent. The inset shows the radial dis- 
tribution functions p'^^\r). Red lines correspond to Eqs. (10) 
and (14). 



on extensive numerical simulations for d = 1, 2, 3, Ferenc 
and Neda [11] proposed 



p{s) ^ n„(s) 



r(a) 



(14) 



where a = (3(i+l)/2. Note that for d = 1 we recover 
Eq. (5). The agreement between Eq. (14) and the nu- 
merical results for P{s) is excellent, as seen in Fig. 2. 



III. NON-POISSONIAN VORONOI CELLS 

More complicate behavior arises when the centers are 
correlated in some way. Although there are few stud- 
ies [3, 7-9] about such systems, the NPV patterns can 
be used to describe qualitatively and quantitative many 
different systems, as we shall see. 



A. One-dimensional non-Poissonian Voronoi cells 

To generate a one-dimensional NPV set of centers we 
proceed as follows. Let p^{x) be the probability density 
to put a new center inside a gap with a scaled size x. 
In a similar way, p^{x) is the probability density that 
the new center is placed at the position x with respect 
to the center at the left of the gap. This kind of model 
proved fruitful in studying the spatial structure of the 
one-dimensional point-island model for epitaxial growth. 
In fact, a suitable choice of p^{x) and p^{x) leads to an 
excellent description of the physical properties of this sys- 
tem [8, 9]. The gap size distribution, p^^\x)^ was shown 



there to satisfy the equation 

^ c^p(Q)(x) ^2^(0). X ^ _p^(^) + 2p^(x), (15) 
dx 

In particular, the probability to put a new center inside 
a gap has the form 



p^{x) = ^p^^\x), 

^7 



(16) 



with ji^ the moment oip^^\x). The probability den- 
sity p^{x) can be written as 



p^(x)= / dz-p^{z)p^{x/z), 

J X ^ 



(17) 



where p^{x/ z)dx/ z is the conditional probability that, 
given a particular gap with size the new center is 
placed inside [x^x + dx]. Unfortunately, in most cases 
the integro-differential given by Eqs. (15), (16) and (17) 
cannot be solved analytically. However, this system can 
be simulated numerically without major difficulties [20]. 

Nonetheless, some exact results can be extracted from 
Eq. (15). If we choose 7 = 1 and p^{\ = x/ z) = 1, the 
probability to put a new center onto an empty site is the 
same for all empty sites on the lattice. Then, with this 
selection of 7 and p^(A), we recover the one-dimensional 
PV case discussed previously in Sec. II-A. As mentioned, 
the position of the centers in the PV case are totally 
uncorrelated, and p^'^\x) is given by Eq. (2), while P{s) 
is given by Eq. (5). 

In general, for 7 > 1 and p^{x/z) = 1, the solution of 
Eq. (15) can be calculated as follows. The simple form of 
the kernel p^{X) leads to dp^{x)/dx = —x^ 
Then, differentiating Eq. (15) we have 



-Vo)(x)//i, 
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dp^ix) 
dx 



+(2 + 7) p^^^(x) = 0. 



(18) 

After some algebra the above equation can be written as 

dp^'^Xx) 
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whose general solution is 
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with /i^ = [r(l/7)/r(2/7)]^ /7. In this case p^{x) is also 
given by Eq. (20). Note that p(^)(0) ^ 0. For this case 
higher spacing distributions, in particular P{s)^ cannot 
be calculated easily. 

Consider now a more general case where p^{x/z) de- 
pends on X and z. We restrict our work to functions 
which are symmetrical about A = x/z = 1/2. This sym- 
metry property comes from the fact that in the absence 
of an external drift (e.g., a field), p^{l — x/z) must be 
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equal iop^{x/ z). Furthermore, we impose the additional 
condition p^(0) = = 0. This property implies that 

the probability to place a new center near an existing one 
is small. 

For large values of dp^^\x)/dx is negative. Then, 
in this regime the behavior of p^^\x) is dominated by 
p^{x). Consequently, for x ^ 1, Eq. (15) takes the form 



dp^^\x) 



dx 



(21) 



which implies p'^^^x) oc exp (— x'^/(7//^)) and p''^\x) 
p^{x). A first correction to this formula can be obtained 
by using the ansatz p^^\x) cx /(x) exp (— 0:^/(7 /i^)) in 
Eq. (15). This procedure gives the differential equation 



(22) 



We conclude that p^^\x) cx x ^exp (—0:^/(7 /i^)). Then 
the behavior of p^^^ (x) for large values of x is completely 
determined by the parameter 7. Furthermore, in the 
limit X 1 the solution of Eq. (15) does not depend on 
p^{X). However, p^{X) controls the behavior of p^^\x) 
for small values of x. In general, for the kind of func- 
tions considered here, for A <C 1 we have p^{X) ~ A^, 
with ( a constant which depends on the functional form 
of p^(A). A series expansion of Eqs. (15) and (17) shows 
that p^^\x) ~ x^ and p^{x) ^ x^ . It is clear that the 
vanishing condition imposed on p^{X) for A = leads 
to p^^\0) = 0. The effective entropic "repulsion force" 
between centers is determined by the parameter given 
by the series expansion of p^{X) around A = 0. 



B. Two-dimensional non-Poissonian Voronoi cells 

In the case d = 2 we proceed as follows. Let q^{s) 
be the probability density to put a new center within a 
Voronoi cell having a scaled area s. Explicitly, we con- 
sider the general form 



on the shape of the Voronoi cell. For example, in the case 
of a circular Voronoi cell with scaled area s, we have 



S + 2 s s 
-TT^r . 



(25) 



The simplest case is 7 = 1 and (5 = 0, for which q'^{r,s) oc 
1/s. In this case every empty point of the lattice has the 
same probability to receive a new center. This, of course, 
corresponds to the PV case discussed in Sec. II-B and 
shown in Fig. 2. 

From a Taylor expansion of Eq. (6) around = 0, 
it is clear that the behavior of p^^\R) is controlled by 
the functional form of c{R). For small values of it is 
reasonable to suppose that the concentration of particles 
within a distance R and R + dR from a given center 
is proportional to the product of g^(R, s) and the area 
dA = 2 7rRdR; then 



c{R) dR - q^'i'R, s) dA - R^^^dR. 



(26) 



Then, p^^\R) - and S controls the effective "re- 

pulsion force" between centers. Using this simple result 
in Eq. (11) we conclude that, for < 1, ^(i?) i^^. On 
the other hand, the behavior of (R) for large values 
of R depends on the behavior of the integral /o^<^^c(^). 

The equivalent of Eq. (15) for the 2D case cannot be 
written easily in terms of known quantities. Following 
Ref. [21] , the effect of a new center on P{s) can be written 
as 



dP(s) ^ , 
as 



M (s) -Mp'^{s)+p'' {s) , (27) 



where M is the average number of preexisting Voronoi 
cells overlapped by the Voronoi cell generated by the new 
center; ^+(5) is the probability density that the new cen- 
ter reduces the Voronoi cell size of a preexisting cell to 
s, p^(s) is the probability density that the Voronoi cell 
of the new center overlaps a preexisting cell with size s; 
and p*{s) is the probability density that the new Voronoi 
cell has size s. 

From their definitions, q^{s) and p^{s) are clearly 
related: q'^{s) takes into account the destruction of a 



(23) 



where jl^ is the 7*^ moment of P{s). In a similar way, 
we define g^(r,s) as the probability density that, for a 
particular cell with scaled size 5, the new center is lo- 
cated at a position r with respect to the center of the 
preexisting cell. For the sake of simplicity, we consider 
just the isotropic case where (with r = |r|) 



g^(r, s) . 



(24) 



In this simplification, the probability to put a new center 
inside a cell depends on the cell itself, regardless of the 
positions of neighboring centers or the areas of their sur- 
rounding cells. The functional form of g^(r, s) depends 
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FIG. 3. In (a) is the initial configuration of centers, while in 
(b) the effect of a new center on the preexisting Voronoi cells 
is shown. Note that in this particular case, the new Voronoi 
cell is almost completely defined by the first three neighbors 
of the new center. 



5 



Voronoi cell by the direct impact of a new center on 
that cell, while p^{s) is more general and expresses that 
the Voronoi cell of a new center overlaps on average 
M preexisting Voronoi cells. Then we can expect that 
p^{s) ^ q^{s) ~ s^P{s). On the other hand, for large 
values of 5, dP{s)/ds < 0, which means that p^{s) con- 
trols the right side of Eq. (27). Thus, in this limit we 
can expect that P{s) ~ exp (—M 5^/(7 /i^)). In our 2D 
model, the tail of P{s) for large values of s depends on 
the parameter 7. 

In the opposite limit s <C 1, the behavior of P{s) is 
given by because this distribution dominates the 



right side of Eq. (27). Thenj9*(s) 



~ s 



Ci 



implies P{s) 



^ s 



C 



In order to understand the behavior of p*{s) for small 
values of s we proceed as follows. 

In the simplest case a new Voronoi cell is completely 
defined by just three nearby centers; see Fig. 3. Note 
that the new center shown in Fig. 3(b) generates a small 
Voronoi cell even though it was formed at the expense 
of the three large Voronoi cells shown in Fig. 3(a). If Pc 
is the probability to have the initial configuration shown 
in Fig. 3(a), then the probability to place a new center 
as in Fig. 3(b) is given by PcA^ where A is the area of 
the target region. Naturally, A scales with the distance 
between centers as A ~ r^. In the case of noncorre- 
lated centers, Pc ~ p^^\ri)p^^\r2 — ri)p^^\rs — r2) ^ 
c(ri)c(r2 — ri)c(r3 — r2) which leads to Pc ^ r^. In the 
PV case we have 7 = 1; thus, p*{s) ^ s'^-^. Note that 
this result agrees with Eq. (14). The case of correlated 
centers is significantly more complicated because Pc can- 
not be written as an independent product of p'^^^s). In 
any case, it is clear that p*(s) for small values of s is re- 
lated with the concentration of centers c{R) through the 
radial distribution functions which in turn in our model 
depend on 6. Then, the value of ( in P{s) ~ for a 
given value of 7 increases with S. However, a general 
relation between S and ( remains unknown. 

As noted, Eq. (27) cannot be solved analytically; how- 
ever, the numerical simulation of this system can be done 
without major complications [22]. 



IV. APPLICATIONS OF THE NPV PATTERNS 

A. Gap size distribution of parked cars 

Rawal and Rodgers (RR hereafter) measured the size 
of the gaps between adjacent parked cars [23]. The data 
(500 gaps) were gathered from four connected streets in 
London without any side streets or driveways. They 
found that small and large gaps are unlikely. The ef- 
fective repulsion between adjacent cars arises because 
drivers need to leave some space between cars to allow 
exit maneuvers. RR suggest that large gaps are unlikely 
because people try avoid the waste of space. However, we 
believe that it is more reasonable to think that this hap- 
pens because people usually prefer to park in large gaps 
where the entry maneuver is the easiest. This implies 
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FIG. 4. (Color online) NPV case with 7 = 4 and p^{X) = 
2 sin^(7rA). Here empirical data A and B correspond to the 
measurements reported in Refs. [23] and [25], respectively. 



that large gaps are often destroyed before small ones. 
RR developed two different models to describe p^^\x); 
however, just one of them describes the statistical be- 
havior of the system properly. In their improved model 
they consider two different factors: people who park any- 
where and those who perform an additional maneuver to 
avoid the waste of space. These assumptions give a good 
description of the empirical data. 

On the other hand, Abul-Magd [24] used the Wigner 
surmise (WS) with /3 = 2 as an approximate model for 
p^^^(x). The WS describes the gap distribution of a one- 
dimensional Coulomb gas at an inverse temperature (3. 
In this system, particles are free to move around a circle 
but experience a logarithmic potential interaction. He 
selected this value of the inverse temperature /3 because 
in this case the WS describes excellently the statistical 
behavior of chaotic systems without time-reversal sym- 
metry, as is appropriate for the car-parking problem. The 
WS does give a good quantitative description of the car- 
parking problem; however, the physical interpretation of 
the logarithmic interaction potential among cars is not 
entirely clear. In Ref. [25], Seba reported more exten- 
sive (over 1200 gaps) empirical data for the gap sizes. 
He modeled the car parking with a Markov chain where 
the cars are allowed to get in and get out of the gaps 
following prescribed probabilistic rules. 

From our model to generate one-dimensional Voronoi 
cells, we can interpret the car-parking problem more sim- 
ply and intuitively. As in the models mentioned before, 
our NPV model for car parking contains some simplify- 
ing assumptions for an ensemble of drivers such as homo- 
geneity and isotropy. Our goal is to analyze a tractable 
version for the problem rather than contend with all the 
subtleties. As mentioned previously, we start with the as- 
sumption that people prefer to park in large gaps rather 
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than in small ones. This follows from the fact that the 
parallel parking is easier in spots where there is space to 
spare. Thus, it is reasonable to propose that the prob- 
ability to park in a gap with size x can be modeled by 
Eq. (16). Naturally we have to choose the appropriate 
value of 7. Once the gap has been selected, the driver 
has to choose the exact parking place inside the gap. For 
gaps shorter than two car lengths the most likely place 
to park would seem to be the middle of the gap. For lack 
of a better simple ansatz, we extend this assumption to 
gaps of arbitrary length. Additionally, the driver should 
avoid parking too close to the cars on the borders of the 
gap in order to guarantee enough space to leave the gap 
when necessary. This implies p^{0) = p^{l) = 0. Fi- 
nally, we claim that p^{X) = p^{l — A); i.e., the drivers 
do not have preference to park near to the car on either 
side of the gap. A simple function which satisfies those 
properties is p^{\) = 2 sin^(7rA). From our numerical 
results we found this functional form for p^{X) gives a 
reasonable description of the empirical data when 7 = 4. 

As shown in Fig. 4, our model describes excellently the 
empirical data given in Refs. [23, 25]. Of course, more 
refined models can be developed, but the most important 
result is that a suitable selection of 7 and p^{\) leads to 
an excellent description and interpretation of the statis- 
tical behavior of the gap size of parked cars. From our 
previous discussion, it is clear that p^^\x) ~ with 
C = 2 for X <C 1. The value of C is related to the ef- 
fective "repulsion force" between adjacent cars. In the 
limit X > 1, p^^\x) ^ x'^ exp(-x'^/(4//^)). Note that 
7 = 4 represents the strong preference of drivers to park 
in large gaps rather than in small ones. 

In the problem of parked cars, we can interpret the 
Voronoi size distribution P{s) through the next-nearest 
spacing distribution p^^\x) as follows. From its defini- 
tion, p'^^\x) clearly gives the probability density that a 
parked car has a spot of length x to perform the exit 
maneuver. From Eq. (4), the distribution of the Voronoi 
cell sizes s is proportional to the distribution of distances 
X that the drivers have to perform an exit maneuver. 



B. Point-island model for epitaxial growth in ID 
and 2D 

In the point-island model for epitaxial growth, atoms 
are deposited onto a substrate where they perform ran- 
dom walks. In the simplest case, when two atoms meet 
they form a static island. In the same way, the atoms 
which reach an island are captured and remain attached 
to it. This case corresponds to irreversible attachment 
and is usually called = 1" in the literature. Another 
important characteristic of the point-island model is the 
fact that the islands do not grow laterally. The size of 
a particular island is given just by the number of atoms 
which belong to it. This system exhibits a scaling regime 
in the limit = F/Z) ^ 00, where F is the deposition 
rate of atoms and D is their diffusion constant. This 



model is, of course, a simplification of the real system 
but it contains most of the relevant physical properties 
required to describe the processes behind the island for- 
mation in epitaxial growth [8, 9, 21, 26, 29-32]. In fact, 
the widely used point-island model gives very accurate 
results in early stages of growth (low coverages) and it is 
an important theoretical model in our knowledge about 
epitaxial growth. 

In this context, the point islands determine the pattern 
of Vononoi cells, playing the role of the centers defined 
in Sec. I. The atoms inside a particular Voronoi cell 
are usually captured by the center of the cell (island). 
Because of this, the Voronoi cells are called capture zones 
(CZs) in the context of epitaxial growth. Naturally, the 
growth rate of an island is related to the size of its CZ. 
For more details about these kinds of models see, for 
examples, Refs. [8, 9, 26-34]. 

1. One- dimensional case 

Consider now the case of a one-dimensional substrate 
with irreversible attachment. A suitable choice to de- 
scribe the spacing and the CZ distributions of the ID 
point-island model is [8] 

/(A) = 30A'(1-A)' (28) 

and 

Blackman and Mulheran [9] originally calculated 
Eq. (28) by first obtaining the average density of atoms, 
ni(x,7/), inside a gap of length y from its expression in 
the stationary state. From this approximation and as- 
suming that the probability of a new nucleation at x is 
proportional to ni(x,?/)^, they found Eq. (28). On the 
other hand, Eq. (29) is based on the numerical results re- 
ported in Ref. [8]. In Fig. 5(a), the results of this model 
are shown and compared with the direct numerical sim- 
ulation of the island nucleation for three different values 
of 5ft; the agreement is excellent. This selection of ^{x) 
and p^{\) reproduce the statistical behavior of the ID 
point-island model with irreversible attachment. For ad- 
ditional information see Refs. [8, 9]. 

2. 2D in circular- cell approximation 

The point-island model with irreversible attachment in 
two dimensions can be modeled following a similar proce- 
dure to that was used previously to described nucleation 
in one dimension. To make analytic progress, we follow 
Refs. [35-38] by approximating the Voronoi cells as circu- 
lar. (While Figs. 1, 3, etc., show that typical cells are far 
from circular, the approximation is better than might be 
anticipated because the average shape over an ensemble 
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FIG. 5. (Color online) The statistical behavior of irreversible nucleation in ID is shown in (a). The CZ distribution for 2D 
nucleation is shown in (b). In all cases our results are compared with direct numerical simulations of island nucleation for 
several values of ?R:. 
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FIG. 6. (Color online) Behavior of the first radial distribution 
p^'^'*(r), which is related to the concentration of islands c(r). 



of cells does tend toward circular.) As shown in Eq. (B2a) 
of Ref. [35] and Eq. (9) of Ref. [38], the isotropic steady- 
state solution of the appropriate diffusion equation with 
flux F inside such a circular Voronoi cell with radius Rc^ 
for a concentric (non-point) island of radius Risi < Rc^ 
gives the following expression for the density of atoms: 



In 



R 

Risi 



1^ 

2 m 



1 



, (30) 



where R (assuming Risi < R < Rc) is the distance from 
the cell center. Note that p{Risi) = (thence increasing 
linearly with R — Risi initially) and dp{R) / dR\^ = 0. 
The first condition comes from the density of atoms being 
zero along the island boundaries. On the other hand, the 
probability of nucleation is proportional to p^{R)^ which 
is maximal along the boundary of the CZ. Hence, p{R) 



is also maximal at = i?c7 as implied in the second 
(Neumann) condition. 

In this framework [with nucleation ex p^(i?)], the prob- 
ability to have a nucleation event within a particular cir- 
cular capture zone with radius Rc^ Pn{Rs), can be writ- 
ten as 



Pn{S) 



P{S) 



(X 



Risi 



dR2TrRp^{R) oc i?^ 



(31) 



Note that Pn{S) is proportional to the third power of the 
scaled area of the cell, i.e., 7 = 3. This is, of course, a 
strong approximation. It is well known that the radius 
of a capture zone fluctuates significantly around its mean 
value because the centers are usually not at the geometric 
center of the cells. Approximating the shape of a CZ by 
a circle neglects those radial fluctuations. Furthermore, 
the density of atoms inside a CZ depends on the position 
of neighboring centers, which is not taken into account 
by Eq. (30). This also implies that the isotropy assump- 
tion used to write Eq. (24) is poor. From Eq. (31) it is 
clear that the probability of nucleation increases with the 
distance from the center and reaches its maximum along 
the boundary of the capture zone. There are many ways 
to select the place of nucleation inside a particular cap- 
ture zone [36-38]. However, in order to keep our model 
as simple as possible, we assign the same probability to 
all points inside a particular CZ, regardless of their dis- 
tance from the center; i.e., a = 0. While this is a crude 
approximation, we can see in Fig. 5 that it is adequate 
to describe P{s). However, previous simplifications have 
important effects on the radial distributions. Figure 6 
shows the behavior of p^^\r). Not unexpectedly, our 
simplified model does not describe p^^^ (r) appropriately; 
i.e., it is not a good approximation for the island density 
c{R). From Eq. (30) it is clear that the concentration of 
islands in the limit <C 1 is given by c(r) ~ R^^ which 
implies S = 1. 
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In order to improve our model, we must take into ac- 
count the fact that p{R) vanishes along the boundaries of 
the islands; i.e., q^{0^s) = 0. Perhaps the simplest way 
to accomplish this goal for the point-island model is to 
propose 

{i ii-Rf<-R<\ 

where < < 1 is a constant, and Rc = (S'/tt)^/^ is the 
average radius of the CZ. From our numerical experimen- 
tation we estimate k = 0.3. In this way, the nucleation 
probability inside a capture zone grows linearly with R 
for points near the island, while it becomes constant for 
points far away. 

Figures 5(b) and 6 also include the results of this last 
model. The description of P{s) is again excellent, but 
now the nearest- neighbor radial distribution is also well 
fitted. Thus, in our model for the island nucleation in 2D, 
the repulsion between centers given by a seemingly has 
a great impact on p^^"^ (r) but is not crucial to determine 
Pis). 

C. Size distribution of second-level administrative 
divisions 

The polygons formed by county boundaries in the size- 
division model resemble Voronoi cells [6, 39]. Inspired 
by this fact, we explore the possibility to use our NPV 
model to study the formation of second-level administra- 
tive divisions (SLADs), such as counties in the USA or 
(non- urban) districts {arrondissements) in France. The 
formation of SLADs depends on many political, cultural, 
ecological, and geographical factors [6, 40]. These fac- 
tors are in general difficult to include in a mathematical 
model. However, our model allows us to give a simple 
interpretation of the SLAD formation in terms of q^{s) 
and g^(r, s). 

We focus here on the results reported by Le Caer and 
Delannay (LCD) for the SLADs in France [6]. The de- 
partments {departements) comprise the first administra- 
tive division in France. In mainland France there are 94 
departments. Each department is divided into several 
districts; there are a total of 322 districts in mainland 
France. Each district has a chief town, mostly with the 
same name as the district. 

LCD compared the area distribution of the French 
SLADs [6, 42] to the one-parameter gamma distribution 
^a{s) with a ~ 4.4 [43]. As seen in Fig. 7(a), the agree- 
ment between 114.4(5) and the data is satisfactory. Note 
that a = 4:A is larger than the value used to describe 
the size distribution of PV cells {a = 2.5). In our model, 
this implies S > 0. It follows that the SLAD formation 
process in the rectangle W was not a completely random 
process. In fact, taking S = 0.5 and 7 = 1 we found 
good agreement between the data and the results of our 
model. The value of 7 was chosen to have an exponential 
tail for P{s) such as in Eq. (14). 



In particular LCD focused on the pattern formed by 
the chief towns of the districts which are located within 
a 544 km x 636 km rectangle whose corners are close 
to the towns Saint-L6, Thionville, Apt, and Mont-de- 
Marsan. The average distance between the enclosed 188 
chief towns is 29.6 km. As shown in Fig. 7(b), the pattern 
formed by the chief towns looks similar to that in Fig. 1. 
However, as noted by LCD, the formation of these SLADs 
was not a PV process. 

The parameters S = O.b and 7 = 1 do not describe the 
cumulative nearest- neighbor distribution, 

F{r) = f dip^^\il (33) 

nor the pair correlation function, 

5W = ;^Er^,Wr + r.-r.)) (34) 

as seen in Fig. 8. Increasing the value of 5 improves 
the estimation of F{r) and g{r) but it leads to a poor 
description of P{s). This means that a model like ours 
is insufficient for this kind of system. 

We attribute this discrepancy to the assumption of 
point islands. One must account for the actual areas 
of cities, which produces an effective short-range repul- 
sion. The simplest way to incorporate this feature into 
our model is to introduce an excluded area around each 
city center in the form of a hard-core radius Vcore- In 
particular, we modify Eqs. (23) and (24) as follows: 

^c(,) ^ (^-Y'°'"^)' p(.)e(. - TrrLe), (35) 

and 

g'^(r, S) (r - Vcoref ©(^ - ^core), (36) 

where 9(^) is the unit step function. When Vcore = 0, we 
recover the original model. From numerical simulation of 
this improved model, we found excellent agreement with 
the data by using Vcorel (^) ^ 0.4, ^ = 2, and 7 = 2 
(see Figs. 7 and 8). This improved model still describes 
adequately P(s), but now the fits involving p^^^{r) and 
^(r) are significantly better. Note that this model does 
not describe the cluster given by Paris and the nearby 
towns. Clearly, the population density is highest in Paris 
and its surroundings. Consequently, it is reasonable to 
expect the formation of clusters of towns in this region 
of France, which is inconsistent with our ansatz of an 
excluded area. 

Our analysis has been applied to SLADs in some 
20 countries, with results generally consistent with the 
above analysis [41]. As we report in Ref. [41], however, 
there are some subtleties and a rich range of nuances, 
e.g., regional differences in the distributions of the areas 
of counties in the USA. 
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(a) (&) 



FIG. 7. (Color online) (a) Size distribution of the 322 districts of mainland France (excluding Paris), (b) NPV pattern generated 
by the 188 chief towns in the districts in a rectangle W in France, as discussed in Ref. [6]. The points A, B, C, and D represent 
Saint-L6, Thionville, Mont-de-Marsan, and Apt, respectively. (The cluster of points in the upper part of the figure represents 
Paris and the nearby chief towns.) The area distribution of these Voronoi cells is included in panel (a) and is rather similar to 
the distribution of actual districts for 1/2 < s < 2. The fits are done iteratively. The improved NPV model [cf. Eqs. (35) and 
(36)] takes into account the finite area of the chief towns by assuming a core radius that is 2/5 of the mean radius. 
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FIG. 8. (Color online) (a) Cumulative distribution F{r) and (b) pair correlation function g{r), based on the set of points 
depicted in Fig. 7(b), using the point-island NPV model and the improved NPV model incorporating the finite size of chief 
towns. 



D. Capture zones of Paris Metro station 

In analogy with the island nucleation, we can define the 
Voronoi cells or CZ for Metro stations as follows. Each 
Metro station represents a center. The Metro stations 
are in competition for passengers in the same way that 
the islands compete for atoms. If we suppose that all 
Metro stations are accessible from anywhere, then most 
of the passengers within a particular Voronoi cell will be 
"captured" by the center of this cell. Of course, this is 



an oversimplification. The passengers not only selected 
a Metro station because of its proximity. They also take 
into account ease of access to it (parking, bus routes, 
commuting possibilities, road conditions, etc.). 

As with the SLADs, in order to have a good set of data 
to apply our model, we have to select a city with a near- 
uniform geographical profile and with a large number of 
Metro stations. These conditions are approximately sat- 
isfied by the capital of France, Paris. Figure 9(a) shows 
a scale diagram of the Metro stations of Paris (RER 
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FIG. 9. (Color online) (a) Voronoi cells (CZ) for the Paris Metro stations. The complete network of stations on the left is 
enlarged to show the stations and associated cells near the heart of the city, (b) Capture-zone distribution of the Paris Metro 
stations. 



stops are also included). Clearly the density of centers 
(Metro stations) is not constant. Naturally, there are 
more Metro stations near the center of the city, where 
the population density is highest. Because of this, the 
largest Voronoi cells are near the outskirts of the city. 
However, we can expect that our NPV model works well 
if we just consider stations near the city center, where the 
density of Metro stations does not change dramatically, 
as seen in the enlargement in the figure. 

For economic reasons it is unlikely that two or more 
Metro stations are very close together. We then expect 
a > because there is an effective "repulsion force" be- 
tween stations. In Fig. 9(b), we compare the empirical 
CZ distribution of the in-town Metro system with our 
NPV model. Good agreement is found with 7 = 2 and 
a = 1.5; the gamma function Ils{s) is also included in 
the comparison presented in Fig. 9(b). Apparently the 
"repulsion force" between Metro stations is bigger than 
that for island nucleation in 2D. 



V. CONCLUSIONS 

Our proposed NPV model can be used to describe and 
interpret many different systems in terms of two inde- 
pendent distribution functions. In ID, p^{x) gives the 
probability density to put a new center inside a gap 
with a scaled size x. This distribution is given explic- 
itly by Eq. (16) and is controlled by the parameter 7. 
This parameter determines the behavior of the gap size 
distribution for large values of x. In fact, for x 1, 
p^^\x) ~ ylx~^exp (— x^/(7/i^)), independent of the 
kernel p^{X). Additionally, 7 modulates the size depen- 
dence of the destruction of gaps. For 7 = 0, the proba- 
bility to put a new center within a gap is the same for all 
gaps, regardless of their size. The larger 7 is, the greater 
is the probability of destruction of large gaps. For the 



car-parking problem we use 7 = 4, which reflects the 
preference of drivers to park in large gaps rather than 
small ones. For island nucleation in ID it is necessary to 
take into account that 7 is a function of the scaled size 
s. 

In ID, the behavior of p^^^(s) in the limit s <C 1 
is completely determined by the fragmentation kernel 
p^(A). For the kernels considered here, we always have 
the generic behavior p^{X) ~ A*^ for A <C 1. The pa- 
rameter ( controls the effective repulsion force between 
centers. For island nucleation in ID and the car-parking 
problem, we found C = 2. For the first system the value 
of ( is fully determined by the density of atoms inside 
the gap in the aggregation regime. In the car-parking 
problem, ( reflects the need of the drivers of allow space 
between cars to perform an exit maneuver. 

In 2D, the probability density q^{s) to put a new cen- 
ter inside a Voronoi cell is also controlled by the param- 
eter 7. In the case of the SLADs, 7 = 1 gives a good 
fit of the empirical data. Because of this, the P{s) of 
many SLADs can be approximated by using the single- 
parameter gamma distribution 110,(5) [40]. In the case of 
Paris Metro stations we used 7 = 2. 

The position r of a new center inside a particular 
Voronoi cell in 2D is determined in our model through 
P^{t^s). For the sake of simplicity, we consider only 
isotropic cases. This probability is closely related to the 
concentration of centers c{R). For example, in the case of 
island nucleation in two dimensions, a good description 
of p'^^^ (r) requires taking into account that the density of 
atoms vanishes along the island boundaries, even though 
the CZ distribution can be well described without taking 
into account this fact. This suggests that many different 
fragmentation models can be used to describe the CZ dis- 
tribution of islands in epitaxial growth. However, just a 
few of them fully describe the statistical behavior of the 
system in a proper way. An additional example is given 
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by the SLADs in France. Two different sets of parame- 
ters describe P{s) properly but just one of them gives a 
good fit for F{r) and g{r). 

We defined the CZs of Metro stations; an extension 
to other systems defined by gas stations, public schools, 
coffeehouses, post offices, etc., is straightforward. In epi- 
taxial growth, the CZ of an island is related to the islands 
rate of capturing atoms. It is reasonable to expect that 
the number of passengers entering a Metro station or the 
influx of customers patronizing a retailer is intimately 
related to the size of its CZ. 
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APPENDIX: RELATION BETWEEN p^^\R) AND 

c{R) 



The model presented here allows us to describe quanti- 
tatively and qualitatively many systems based on simple 
assumptions about them. For example, in the car park- 
ing problem we based our model on some assumptions of 
the driver preferences related with parallel parking. Our 
ansatz leads to a reasonable description of the gap size 
distribution. For the description of the CZ distributions 
in the point-island model, we based our simulation on an 
estimate of the density of atoms (besides other observa- 
tions) , which comes from the direct numerical simulation 
of this system. Nevertheless, there are systems where the 
NPV patterns are determined by factors difficult to es- 
tablish as a mathematical expression. In the case of the 
Paris Metro stations, it is clear that there is an "effective 
repulsion" between stations because it is not economi- 
cally viable to put two or more stations too close. In 
our model this implies ^ > 0. However, there are other 
political, historical, and geographical factors which also 
affect the CZ formation. Something similar happens in 
the case of the SLADs. 



Despite its implicit simplifications (such as homogene- 
ity and isotropy), our model proves to be a powerful tool 
to describe several complex systems which are defined 
through an array of points. 



For an arbitrary isotropic concentration in a circle of 
radius Rq, the expected number of particles inside a con- 
centric disk with radius R < Rq is dr 27Tr c{r) . The 
probability that this disk contains no particles is 



V Jo drrc{r)_ 



fg° dr2wrc(r) 



when i? <C -Rq- This last statement leads to 



/o drrcjr) 
J^drrc{r) 



< 1. 



Then, we have (for i? <IC -Rq) 

E{R) = e-/o''rf'-2'rrc(r)_ 



(Al) 



(A2) 



(A3) 



By definition, E{R) is the probability to have an empty 
disk with radius R; then E{R) — E{R-\-dR) is the proba- 
bility of have an empty disk with some particles between 

and -h dR] i.e., p^^\R)dR. We conclude that 



p^'\R) 



dE{R) 
dR ' 



(A4) 



Equation (6) can be obtained from Eqs. (A3) and (A4). 



[1] D. Weaire, J. P. Kermode, and J. Wejchert, Phil. Mag. 

B 53, LlOl (1986). 
[2] D. Weaire and N. Rivier, Contemp. Phys., 25, 59 (1984). 
[3] P. A. Mulheran and D. A. Robbie, Europhys. Lett., 49, 

617 (2000). 

[4] F. Mercier and O. Baujard, Voronoi diagrams to model 
forest dynamics in French Guiana, in 2^^ International 
Conference on GeoComputation, University of Otago, 
NZ (1997). 

[5] M. Xue, Airspace Sector Redesign Based on Voronoi Di- 
agrams, Proc. AIAA Guidance, Navigation and Control 
Conference, Honolulu, HI (2008). 

[6] G. Le Caer and R. Delannay, J. Phys. I (France) 3, 1777 
(1993). 



[7] L. Zaninetti, Phys. Lett. A, 373, 3223 (2009). 

[8] D. L. Gonzalez, A. Pimpinelli, and T. L. Einstein, Phys. 

Rev. E 84, 011601 (2011). 
[9] J. A. Blackman and P. A. Mulheran, Phys. Rev. B. 54, 
11681 (1996). 
[10] T. Kiang, Z. Astrophys., 64, 433 (1966). 
[11] J.-S. Ferenc and Z. Neda, Physica A 385, 518 (2007). 
[12] H. J. Hilhorst, J. Stat. Mech. P09005 (2005). 
[13] H. J. Hilhorst, J. Stat. Mech. P05007 (2009). 
[14] T. Xua and M. Lia, Phil. Mag. 89, 349 (2009). 
[15] E. N. Gilbert, Ann. Math. Stat. 33, 958 (1962). 
[16] D. L. Gonzalez and G. Tellez, Phys. Rev. E 76, 011126 
(2007). 

[17] J. P. Hansen and I. R. McDonald, Theory of Simple Liq- 



12 



uids, 2nd ed. (Academic Press, London, 1986). 
[18] D. L. Gonzalez, A. Pimpinelli and T. L. Einstein, e-print 

arXiv:1111.5212 (submitted to PRE). 
[19] S. Redner and D. ben-Avraham, J. Phys. A 23, L1169 

(1990). 

[20] A lattice with N sites is taken and an array of points 
on the lattice is generated. Each new point is introduced 
via weighted random number generation. First, a pair of 
existing points is selected between which to introduce the 
new one. The selection is weighted by the 7*^ power of 
the separation of the points. Then the actual position in 
the gap for insertion, x, is selected with weighting given 
by the function p^{x/y). We use L = 10^ and gather data 
for 100 particles over 5 x 10^ realizations. 

[21] M. Li, Y. Han, and J. W. Evans, Phys. Rev. Lett. 104, 
149601 (2010). 

[22] In the 2D case we use a square L x L lattice with L = 
200. To place each new center, a Voronoi cell is chosen 
weighted by the 7*^ power of the size of the cell. Then 
the position r of the new center is selected by using the 
distribution q^{r,s). To achieve good statistics we use 
p = 5 X 10~^ and gather data over 5 x 10^ realizations. 

[23] S. Rawal and G.J. Rodgers, Physica A 346, 621 (2005). 

[24] A.Y. Abul-Magd, Physica A 368, 536 (2006). 

[25] P. Seba, J. Phys. A 41, 122003 (2008). 

[26] F. Shi, Y. Shim, and J. G. Amar, Phys. Rev. E 79, 011602 
(2009). 

[27] J. W. Evans and M. C. Bartelt, Phys. Rev. B 63, 235408 
(2001). 

[28] M. N. Popescu, J. G. Amar, and F. Family, Phys. Rev. 



B 58, 1613 (1998). 
[29] J. G. Amar and M. N. Popescu, Phys. Rev. B 69, 033401 

(2004) . 

[30] C. Ratsch, Y. Landa, and R. Vardavas, Surface Sci. 578, 
196 (2005). 

[31] V. L Tokar and H. Dreysse, Phys. Rev. B 80, 161403(R) 
(2009). 

[32] T. J. Oliveira and F. D. A. Aarao Reis, Phys. Rev. B 83, 

201405(R) (2011). 
[33] J. W. Evans, P. A. Thiel, and M. C. Bartelt, Surface Sci. 

Rept. 61, 1 (2006). 
[34] P. A. Mulheran, in Metallic Nanoparticles, edited by J. 

A. Blackman (Elsevier, Amsterdam, 2008), Chap. 4. 
[35] J. W. Evans and M. C. Bartelt, Phys. Rev. B 66, 235410 

(2002). 

[36] M. Li and J.W. Evans, Multiscale Model. Simul. 3, 629 

(2005) . 

[37] M. Li, M.C. Bartelt and J.W. Evans. Phys. Rev. B 68, 
121401(R) (2003). 

[38] M. Li and J.W. Evans, Surface Sci. 546, 127 (2003). 

[39] See [http://historical-county.newberry.org/website]. 

[40] R. Sathiyanarayanan, Ph.D. thesis. University of Mary- 
land, College Park, 2009. 

[41] R. Sathiyanarayanan et al., preprint. 

[42] See [http://eii.wikipedia.org/wiki/List_of_ 

arrondissemeiits_of _France] . 

[43] The gamma distribution was also used in Refs. [40, 41] 
to describe the area distribution for the original Thirteen 
Colonies that evolved into the USA. 



